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INTRODUCTION 

Engineering codes are needed which can calculate convective heating rates accurately 
and expeditiously on the surfaces of high-speed vehicles. One code which has P^e™. 
these needs is the Langley Approximate Three-Dimensional Convective Heating (LATCH) cod 
ref n it ui the axisymmeinc analogue in an integral boundary-layer method to calculate 
laminar and turbulent healing rates along inviscid surface streamlines. It requtres the soluuon of 
!“scid now field to provide the surface properties needed to calculate the stream hues and 
streamline metrics The LATCH code has been used with inviscid codes which calculated 
flow field on structured grids. Several more recent inviscid codes calculate flow field properties 
on unstructured grids. The present research develops a method to calculate inviscid surface 
streamlines, the ftreamline metrics, and heating rates using the properties cal ^ ulat ^ d f '°^ ASA 
Sd flow fields on unstructured grids. Mr Chns Riley 

LaRC developed a preliminary code in the C language, called UNLA , P 

these goals No publication was made on his research. The present research extends a 
improves on the code developed by Riley. Particular* mention at, ° 

region, and the method is intended for programming in the FORTRAN 90 lan 0 ua 0 e. 

TNVISCID SURFACE STREAMLINES AND METRICS 


Let R be the position vector 


R = xi + y j + zk 


( 1 ) 


( 2 ) 


The body surface is described by 

F(x, y, z) = 0 

Therefore the Cartesian coordinates x, y. and z are not independent. We will consider )' and z as 
independent and that * can be determmed from eq. (2) as some funchon of y and z. The method 
developed here is similar to that in ref. 2 for the metric. 

The unit normal to the body surface is given by 


VF 


e„ = 


F x i + F y ] + F z i 

2 il/2 


(3) 


" l VF l [Fx + Fy + F z ] 


The velocity is given by 


V = u i + v j + w k 


( 4 ) 


On the body surface, V and e n are perpendicular. Therefore, 

V ° e n =0 
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and hence 


i iF x + vF y + wF z = 0 


or 


u-- 


( 1 vFy + wF z ) 

Tx 


(5) 


We will also use streamline coordinates t and (3, where S = \Vdt is distance along a streamline 
and P is a coordinate constant along a streamline. The differential arc lengths are given by 

ds = Vdt 

/ip dp 


and 

where hp is the metric 
is e s , where 


( 6 ) 

(7) 


along a t = constant line. The unit vector in the direction of a streamline 


V =Ve< 


( 8 ) 


and the unit vector perpendicular to e s and e n is 

ej_ = e s x e n 


(9) 


is therefore tangent to the surface and normal to e s . If we consider t and P as independent 


variables for surface streamlines, then 


dR 

dt 


= V 


and 


dR 

ap 






( 10 ) 


(ii) 


the unit vector tangent to a t = constant line on the surface and it is not 


ep is 

perpendicular to e s . 

We set 

ep = cos r) £j_ + sin tj e s 

where ri is to be determined. Let n be the straight line coordinate normal to the body. Then 


( 12 ) 


4 


dR „ 
dn " 


( 13 ) 


and R = K(r, p, /i) . Also = x l + y j + y k and the derivative of R gives 

dR = dxl + dyj + dzk = Vdte s + /ipdpep + dne n (14) 

We note that t, P, and n are not orthogonal coordinates. Substitute eq. (12) for <?p in eq. (14) to 
get 

dxi + dyj + dzk = Vdte s + /ipdp(cos 
Since y and z are functions of t, P, and n, we can write 


and 


dz , dz 
— dt + — 
dt ap 


ej_ + sin r| e s ) + dne n 

(15) 


(16) 

dn 


— dn 

(17) 

dn 



where t, P, and n are independent coordinates and n = 0 on the body surface. Substitute eq. (16) 
into eq. (15) and take the dot product of both sides of the resulting equation with j . Equate the 
coefficient of dp on both sides to get 


3y 

ap 


= /ip[cos r\ e_L 0 j + sin ri e s ° j] 


(18) 


Similarity, substitute eq. (17) into eq. (15), take the dot product of both sides with k , and equate 
the coefficients of dp to get 


= /ip[cos ri e± ° k + sin r\ e s ° k] 

Since e ± =e s xe n , then use eqs. (3) and (4) to perform the cross product. The result is 

[i(vF z - wFy ) + j(wF x - uF z ) + k{nF y - vF x )} 

= V | VE 


(19) 


( 20 ) 


Using eq. (20) in eqs. (18) and (19) we get 
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and 


dy 

S 

dz 

ap 


COS T| 


V | VF | 


+ sin 


/j P 


cos r| 


(uF, - vF x ) 
V | VF | 


+ sin 



( 21 ) 

( 22 ) 


We can eliminate sin r\ from eqs. (21) and (22) to get 


/ip cos r) = 


dy dz 

w — — v — - 

| VF | 

ap ap^ 



VF X 


(23) 


In obtaining eq. (23), u was found from eq. (5). The left side of eq. (23) is the metric 
h = h$ cos r\ needed for the axisymmetric analogue to replace the axisymmetnc body radius. 

Note that hd(3 is the perpendicular distance between adjacent streamlines as shown in Fig. 1. 



We can also eliminate cos r\ from eqs. (21) and (22) to get 


dy 


dz 


('’F x -uF,,)-y- + ^F x -uF,) — 


/ip sin r| 




ap 


(24) 


F X V 



6 


and again u can be obtained from eq. (5). Equations (23) and (24) allow us to determine ri if 
needed. Equation (23), however, is the relation needed to obtain the metric h = /ip cos r| for the 

axisymmetric analogue. An alternate method is given in the appendix which uses orthogonal 
coordinates. 


Stagnation Region 

We will consider first the stagnation region for a body which is symmetric about the x-y 
plane. The stagnation point is located at y = yq and Z = °- In the re g ion near the stagnation 
point we can approximate v and vv by the relations. 


and 


v ~ v y (y - y 0 ) 


(25) 


VV = W z z 


(26) 


where 


Since 


v v = 


'dv' 
K dy JO 


and vv T = 


( dvv^ 

[dz J ( 


dy , dz 

— = v and — = vv , 

dt dt 


eqs. (25) and (26) can be written as 


and 


dy 

dt 


= v Y (y-y 0 ) 


dz 

ar 


= w z z 


(27) 

(28) 


Equations (27) and (28) can be integrated, holding (3 constant, to give 

y -, 0 =C,(|5)/-' P9) 

and 

z = C 2 ($)e Wzt < 30) 

where C t (P) and C 2 (P) are functions of (3 and thus constant along a streamline. The stagnation 
point y = y 0 and z = 0 requires t = - ~ since v y > 0 and w y > 0 in general. We have some 

liberty in defining Cj(P) and C 2 (P) ; therefore, we set 


Cj ((3) = 8 cos (3 


(31) 
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and 


C 2 (P) = e 



(32) 


V ♦ 

where e is a small parameter with the dimension of length. The factor l-j- in eq. (32) is 

inserted for convenience to make the curve t = 0 perpendicular to the streamlines which cross 
that line. However, all the other lines of t = constant are not perpendicular to the streamlines. 

Using eqs. (31) and (32) in (29) and (30) give 

y - y 0 = £ cos P e Vyt ^3) 


and 



v v 

— sin P e 
vv. 


W z t 


(34) 


Eliminating t from eqs. (33) and (34) gives another form of the streamline equation as 





(v \ 




vv. 




\ " / 

y — yo = £ cos P 


Z 



E, 

V v 

1— sin P 



V 

1 



(35) 


where, again, p is constant along a streamline. If we eliminate p from eqs. (33) and (34) we get 


y ~ >’0 


8 e 


vj 


+ 


- 1 ? 


v y w t 
— — e 4 
vv. 


= 1 


(36) 


which is the equation for an ellipse for t = constant lines. Equations (35) and (36) show that 
streamlines (P = constant) are perpendicular to t = constant lines only for t = 0, unless v y - w z . 

For eq. (23) we can use eqs. (25), (26), (33) and (34) to determine the factor 




r 

1 

— 

dy dz 
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- vv, e , 

t sin 2 P-v e 2 1 
w, V 

^cos 2 P 

ap apj 



W z 
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n v ~ 

Z 2 J v y w z 

sin" P + — cos" P 
w z 


(v v + w.)t 

e y 


(37) 


Then eq. (23) becomes 


VF 




J v y w z 

fi 4 — C0S~ 6 

V 

J l 1 1 p 1 vVJ p 

W z 


(V +W )1 

e - 


(38) 


Note that F x < 0 for the nose region of a blunt nosed body in order to make e n the unit outer 
normal to the body. By using eq. (5) for u and eqs. (25). (26), (33) and (34) the magnitude of the 
velocity becomes 


= y[u 2 + v 2 + vv 2 = 


(F y v + F,w) z 2 2 

— I + v z + w z 


1/2 


= [(F 2 + F y )v~ + (F 2 + F 2 ) w 2 + 2F y F z vw\ 

F x I 


1/2 


F x \o 


(F 2 + F 2 )o cos 2 P e 2V + {Fx + F z 2 )o sin 2 P e 2w *‘ 


VV T 


V v • o O (V y +W,)f 

+ 2(F V F Z ) 0 J— sinPcosPe y 


VV T 


1/2 


(39) 


For bodies symmetric about the x-y plane, F z = 0 in the plane of symmetry. 
For the laminar heating rate, we use dS = V dt to determine the integral 


I = )h 2 Vds = \h 2 V 2 dt 
o -°° 


(40) 


Using eq. (38), eq. (40) can be integreted to yield 


I = 


l Vf lo e 4 v v w. 

sin 2 p + — cos 2 P 

(F x )l ' 

J 


2 ^2(v v +vv,)f 
2(v y + vv, ) 


( 41 ) 


For the laminar momentum thickness, we need the ratio I [l2 /hV which becomes 
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hV ^2(vT + iv.) 


(42) 


Note that this result is independent of (3, t, and e in the stagnation region. 


Region Bevond the Stagnation Region 

It is convenient to let t = 0 at the edge of the stagnation region. Then - °o < t < 0 for the 
stagnation region, where the equations given above are valid. 

For t > 0, we can determine the metric h = cos T] from eq. (23) provided we can 

calculatae — and — . On the body surface y and z are independent variables and also t and (3 

3P dP . . , 

are independent variables. The velocity components are given at specific positions x, y and z 
from the in viscid solver. We note that 


f^ = v(43) 
dt 

d 1 }' _ ( 44 ) 

apa? ap 


But v = v(y, z ) and 



apar 


- L Therefore eq. (44) can be written as 

arap* 


a 

' r)y ' 


r *f| 


' dv ' 

dt 

w 




K d Zj 


(45) 


In a similar manner 


d_(<h' 

dt{d p, 


, dy j. 5P 



dz 

ap 


(46) 


3 3z 

Starting at the ege of the stagnation region, where t = 0 and and are known, we can 


determine y, z, % and || along a streamline by integrating eqs. (43), (45), (46) and ^ = w 
ap ap 

numerically with t or S = J Vdt as the independent variable. In order to do that, we must be able 
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dv dv dw 3w 

to interpolate the inviscid solution on the surface to determine v, vv > ^ - ^7 > 3 ^ and ^7 r 
equations given above. 

Surface Fitting Geometry 

Consider the inviscid solution on an unstructured grid. For the body surface, we would 


know 


Xj , Hi , , Vj , Zi , and vv, 


at points on the surface for Cartesian coordinates (x, y, z) with corresponding inviscid velocity 
components («, v, w). The points do not, in general, correspond to any regular pattern on the 
surface. The equation for the body surface is generally not known, but can be represented by eq. 

(2) Then the unit normal, e n , to the surface, although not known, could be determined from eq. 

_ dF 

(3) where ( 1 , ], k ) are unit vectors in the (x, y, z ) directions, respectively. Also, F x - ^ , 
f — and F, = — . As shown earlier, V ° e n =0 implies 


uF x + vFy + wF z = 0 


We will consider 3 adjacent points on the body surface as shown in Fig. 2. 



Figure 2. Adjacent points on body surface 
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We wish to surface fit the conic equation 


F(x , y,z) = y 2 + &[y + “ii 1 + a 3 z + a 4 x2 + a 5 x + a 6 - 0 

to the region inside the triangle formed by the 3 points. The coefficient of y is chosen to be 
unity, which means this term cannot be zero. We could have chosen some other term to have 
unity as its coefficient, as long as we know that term is not zero. With x the coordinate in the 

longitudinal direction, then the y 2 term would not be zero except for some flat surfaces. 


For the surface fit to the triangular region, we need to determine the coefficients a l ,a 2 , 
a 3 , a 4 , a 5 , and a 6 . Therefore we need 6 equations. These equations follow by applying 
equtions (5) and (47) to the points 1, 2, and 3 on the surface. From eq. (47) we get 


F x = 2 a^x + a$ 
F y = 2y + flj • 
F z - 2a 2 z + a 3 


(48) 


then eq. (5) becomes 


u(2a^x + a$) + v( 2 y + flj) + w(2a 2 z + a 3 ) — 0 


(49) 


Rearranging, we get 


ayv + a 2 (2wz) + a 2 w + < 34 ( 210 :) + a$u - - 2vy 


(50) 


Now apply eqs. (47) and (50) at the 3 points on the surface. This gives 6 equations for the 6 
unknowns a ^, ..., a 6 . After determining these 6 coefficients, we can calculate the surface 

normal, e n , from eq. (3) and .v as a function of y and z from eq. (47). However, there is no 

guarantee that the surface equation is continuous with adjacent triangles to the 3 points. Note 
that F(x, y,z) = 0 being the body surface, there are only two independent variables. It is 
convenient to use y and z as these two indepdendent variables since the streamines near a blunt 

1/2 - 


for „y- 


nose (x = 0 , y = 0 , z = 0 ) has the distance along the streamline varying as y and z but x 
direction. This can be seen more clearly when looking at the stagnation region. As an example, 
consider a blunt-nosed body described by the simplier conic shape 


y 2 + Kz 2 = 2ax + Bx 2 ( 51 ) 

where the cross section for x — constant is an ellipse and K, a, and B are constants. K determines 
the ellipticity for the cross sectional shape, “a” is the nose radius of curvture in the (v - x) plane, 
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(a/K) is the nose radius of curvature in the (z - *) plane, and B determines the conic section in the 
jc-direction, i.e., 

B < 0 for elliptic, B = 0 for parabolic and B > 0 for hyperbolic longitudinal shapes. 


Here we write eq. (5 1) as 


F( x, y, z) = y 2 + kz 2 - 2 ax - Bx 2 = 0 


(52) 


We can solve for jc from the quadratic formula, with B * 0, using x - 0 at y - 0 and z 0, to 
determine the proper sign. Thus eq. (52) becomes 


x = 


-a + [q 2 +g(y 2 +fo 2 )] 1/2 

B 


(53) 


which also can be written as 


+ Bx = [a 2 + B(y 2 + kz 2 )] 


2v,l/2 


(54) 


Equation (54) is convenient because eq. (52) gives 

F = — = - 2a - 2Bx = - 2(a + Bx) = - 2[a 2 + B(y 2 + kz 2 )] 1/2 
* dx 


Also 


7\F 

F v = — = 2y and F z = 2 kz 
y dy 


(55) 

(56) 


For zero angle of attack, the nose point is the stagnation point, and for the region near the 
stagnation point 

v — v v v and w ~ w z z 


(57) 


where v y = 


f 

— and w 

i Jo 

u is determined from eq. (5) as 


l dz J 


u = 


— (vF v + wF z ) 

F y 

1 X 


and using eqs. (55), (56), and (57) we get 
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[v y 2 + w z Kz 2 ] 


u ~ 


[a + B(y~ + Kz )] 


2 mL/2 


(58) 


and near the stagnation point, 


u ~ 



+ vvv 


.Ki- 


el 


(59) 


Thus u ~ y 2 and z 2 whereas v ~ y and vv ~ z. Therefore, we will use y and z as independent 
variables on the surface and * can be determined from F(x , y, z ) = 0. Also, we will work with v 
and w since u can be determined from eq. (5). 


Surface Fitting Velocity Components 

If we only use the 3 points on the triangle 1, 2, 3, then the velocity components v and vv 
can only be fit as linear variations in y an z, i e., 


v = b[ y + ^ + t >2 

w = qy + c 2 z + 


(60) 


where the constants can be determined by applying the equations above for v and vv to the 3 
points where v and vv are known. If we wish more than linear variations in y and z, then we must 

consider more than 3 points. 

One possibility is to add 3 adjacent points as shown in Fig. 3. 



4 


Figure 3. Adding three adjacent points 
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With 6 points, we could write v and vv as 

v = b^y + b 2 z + £>3 + ^4>' 2 + + b^yz 

vv = qy + c 2 z + c 3 + c 4 y 2 + c 5 z“ + c 6 yz (61) 

The 6 coefficients for v and for vv could be determined by applying the equations above to the 6 
points where y, z, v, and vv are known. Although 6 points are used to determine the coefficients, 
the equations above are applied only to the region inside the triangle formed by points 1, 2, and 

3. 

If we go outside the triangle 1, 2, 3 then other points would be used for that region. For 
example, suppose we go beyond (1, 2, 3) into the triangle (2, 5, 3). Then we would add 2 points 
(7 and 8) to the points (1, 2, 3, 5) to make 6 points to determine the coefficients for the region (2, 

5, 3). See Fig. 4. 

\ 

1 

I 

/ 

5 

\ 

\ 

\ 

\ 

-_k 

7 


Figure 4. Adjacent triangle (2, 5, 3) 

What is the continuity of the v and vv along the common boundary line 2-3? Along the line 2-3, 
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y ~ >’2 _ £ ~ z 2 (62) 

^3 ~ yi z 3 ~ z 2 

I 2 

. (y 3 - y2)(z-z 2> 

y 9 H 

(Z3-Z2) J 


Thus V 2_3 = by 


(y 3 - yi) , \ 

> 2 + 7 ^ \(Z~Z 2 ) 

(z 3 “ Z 2 ) 


+ b 2 Z + 63 + 64 


+ ^ 5 z 2 + £> 6 z 


3>2 


+ 


^ ^3 ~ ^2 ^ 
bz-i ~ Z2 


(z - z 2 ) 


(63) 


Let V 2_3 = + ^ 2 Z + #3 for Z 2 - z ^ Z 3 (6 

The coefficients B h B 2 , and 5 3 must be the same for triangle (1, 2, 3) as those for triangle ( 2 , 
5, 3). These conditions would replace pt. 1 for (2, 5, 3). Note that the v 2 _3 equation above 
already satisfies v = v 2 and v = v 3 . Thus the relation given by eq. (63) or (64) along with v at 
points 5, 7, and 8 allow us to determine the 6 coefficients in eq. (61). Similarily for the w 
component of velocity. 

CONCLUDING REMARKS 

The method developed in this report can be used to calculate the heating rates on high 
speed vehicles. The improved solution for the stagnation region should provide more accurte 
heating rates and the technique can be used for regions downstream on unstructured as well as 

structured grids. 
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APPENDIX 

Alternate Surface Coordinates 

The use of t and (3 as independent coordinates on the body surface did not allow t = 
constant lines to be perpendicular to surface streamlines, along which (3 = constant. An alternate 
method is to use 0 and P as independent coordinates where 

dS = Vdt = h s dQ (A-l) 

and 0 = constant lines are perpendicular to the surface streamlines, where p is still constant along 
a streamline. The metric associated with dQ is h s and it is generaly not needed to calculate 

heating rates. The coordinate 0 is defined such that the unit vector in the p direction is e ± . 

Thus eq. (14) becomes 

dR = dxi + dyj + dzk = h s dQe s + hd$ e ± + dne n (A-2) 

Now we are dealing with orthogonal coordinates 0, P, and n, and the metric needed for the 
axisymmetric analogue is h. Considering y and z as functions of 0, P, and n then eq. (A-2) dotted 

with ] gives 


dy , ~ h(wF x -uF) 

ap x v | vf | 


where eq. (20) was used for e ± . Dotting Eq. (A-2) with k gives 


dz ,, r h{uF y -vF x ) 

ap x v | vf | 


(A-3) 


(A-4) 


Note that — is now a partial derivative holding 0 and n constant, whereas the — in eqs. (21) 

ap 

and (22) held t and n constant. Multiply eq. (A-3) by w and subtract eq. (A-4) multiplied by v to 
get 


h = 




VF 


VF r 


(A-5) 


0 

Eq. (A-5) is similar to eq. (23) except here — holds 0 and n constant. 

Consider the stagnation region again, but using 0, P, and n as independent variables. 0 is 
related to t by eq. (A-l) and t must be considered a function of 0 and P on the surface. Eqs. (33), 
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3 

(34), (35), and (36), are valid here also. However, when determining — and ^ from eqs. (33) 
and (34), we get 


(A-6) 


r dy ' 

\ 

= ee 

'y f 

. o ft 9r 

- sin p + v v cos p — 


0 


L y 9P 


and 


f 3z l 

= £ , 

h - e":' 

R ■ R dr 

cos p + w T sm p — 


0 1 

1 

9pJ 


(A-7) 


Interestingly, when we form 


9v 9z 

w— — v — 

9(3 9(3 


dt 

the terms involving — cancel and eqs. (37) and 

9(3 


(38) are also valid here. Although t is not an independent variable, we can still use it and eqs. 
(41) and (42) are also valid. 

For the region beyond the stagnation region, we could use eq. (A-5) to determine h if we 
can determine — and — from adjacent streamlines. On the other hand, if we are tracing only 

9(3 9(3 

one streamline, then we must integrate the differential equation derived below to determine h. It 
follows from eq. (A-2) that 


9 R , „ 

— = n,. 

90 5 * 


(A-8) 


and 


dR , „ 
— = he i 
9(3 1 


(A-9) 


We note that the equality of cross derivatives gives 

9 2 R _ d 2 R 

9(390 ~ 909(3 


(A- 10) 


Differentiate eq. (A-8) with respect to (3 and eq. (A-9) with respect to 0 and then substitute into 
eq. (A- 10) to get 


9e, 9 h 

J i 


+ 


9p 9p 


de ± 9/i . 

e . = h—± + — e I 

5 90 90 ± 


(A-ll) 


Take the dot product of eq. (A-ll) with e L and divide both sides by h s h to get 
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1 dh 1 de s 

= — o e | 

h h s dB h 3(3 

The partial derivative on the right side can be evaluted by using the transformation 


(A- 12) 


1 de, 


f \ 


h ap 


de 


\ ( h' i 

V J Jz,n 


1 dy f de s ^ 
h 3p ^ dz j 


1 dz 


(A-13) 


Eqs. (A-3) and (A-4) give relations to determine — and — . Then eq. (A- 12) can be 

h dp h dp 

integrated along a streamline using 

dS = h s dB 

. h t dB , ...... 

as the variable of integration or with dt = — — as the variable of integration. 


If we dot eq. (A- 11) with e s we get 


3 h s de± 

— - = h — — ° e , 

ap ae 5 


(A-14) 


de, ~ ~ de s 

Divide both sides by h,h and use the vector relationship — — - ° e s = - e ± ° to write eq. 

4 39 


(A-14) as 


i a/i. 


1 de ( 


e± 


39 


h s A dp h s 39 
The partial derivative on the right side can be transformed using 


_i_a^_ 

h s 39 


. 3v . 

V J Jz,n 


1 3y + f de s ) 1 3z 


h r 39 


dz 


V ^ Jy,n $ 


h r 39 


(A- 15) 


(A-16) 


and from eq (A-2) we obtain 


1 3y ^ v 

hr 39 5 V 


(A-17) 



19 


and 


h s 50 


e s ok 


vv 

V 


(A- 18) 


Eq. (A- 15) gives the curvature term 

1 1 dh s 
h s h 5(5 

which is needed when solving the cross-flow boundary layer equation. For heating rates, 
however, the cross-flow boundary layer equation is not needed. 


